%% This Function compute the vector of probability inequality constraints theta in the paper 
function [mu,mu1,mu2,mu3,mu4]=incosp(Y,Z,D,A)
%Estimate the 2 conditional probabilities P(D=1|Z=z)
P11=sum(Z==1&D==1)/sum(Z==1);
P10=sum(Z==0&D==1)/sum(Z==0);

% 
y1=Y(Z==1);
y0=Y(Z==0);
d1=D(Z==1);
d0=D(Z==0);

%Estimate the moment inequality constraints
mu1=zeros(max(size(A))-1,1);
mu2=zeros(max(size(A))-1,1);
mu3=zeros(max(size(A))-1,1);
mu4=zeros(max(size(A))-1,1);

for i=1:max(size(A))-1
    mu1(i)=(sum(y1>=A(i)&y1<A(i+1)&d1==1)/length(y1))-(sum(y0>=A(i)&y0<A(i+1)&d0==1)/length(y0));
    mu2(i)=(sum(y0>=A(i)&y0<A(i+1)&d0==1)/length(y0))-(sum(y1>=A(i)&y1<A(i+1)&d1==1)/length(y1)-(P11-P10));   
    mu3(i)=(sum(y0>=A(i)&y0<A(i+1)&d0==0)/length(y0))-(sum(y1>=A(i)&y1<A(i+1)&d1==0)/length(y1));
    mu4(i)=(sum(y1>=A(i)&y1<A(i+1)&d1==0)/length(y1))-(sum(y0>=A(i)&y0<A(i+1)&d0==0)/length(y0)-(P11-P10));

end


 mu=-[mu1;mu2;mu3;mu4];